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Abstract. This article is concerned with the question of constructing efficient multigrid precon- 
ditioners for the linear systems arising when applying semismooth Newton methods to large-scale 
linear-quadratic optimization problems constrained by smoothing operators with box-constraints on 
the controls. It is shown that, for certain discretizations of the optimization problem, the linear 
systems to be solved at each semismooth Newton iteration reduce to inverting principal minors of 
the Hessian of the associated unconstrained problem. As in the case when box-constraints on the 
controls are absent, the multigrid prcconditioner introduced here is shown to increase in quality as 
the mesh-size decreases, resulting in a number of iterations that decreases with mesh-size. However, 
unlike the unconstrained case, the spectral distance between the preconditioners and the Hessian is 
shown to be of suboptimal order in general. 
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1. Introduction. The objective of this article is to develop efficient multigrid 
preconditioners for the linear systems arising in the solution process of large-scale 
PDE-constrained optimization problems using semismooth Newton methods. The 
model problems of interest have the form 

{minimize ^\\ICu — yd\\^ + f /3 > fixed, 

subj. to : a ^ u ^ b , i^-^) 

where /C : W — >■ 3^ is a linear compact operator between the two function spaces hi and 
3^, and a,b £U,yd & y are given functions. We regard K. as the solution operator of 
a linear partial differential equation (PDE): ii e : y x U ^ y* defines the linear PDE 

e(y,w) = (1.2) 

then e(jj,u) = if and only if y = ICu. This way we obtain an equivalent formulation 
of (|l.ip that has become standard in the PDE-constrained optimization literature [T3] : 

I minimize ^\\y - yaf + §\\uf , 

[ subj. to : e{y, u) = 0, we Z^ad — {u : a ^ u ^ b} . 

In this formulation u is the control and we shall call y the state. We shall also refer 
to p.ip as the reduced form of (|1.3p . 

Due to the availability of increasingly powerful parallel computers, the scientific 
community has shown a growing interest over the last decade in developing scalable 
solvers for large-scale optimization problems with PDE constraints. Multigrid meth- 
ods have long been associated with large-scale linear systems, the paradigm being that 
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the solution process can be significantly accelerated by using multiple resolutions of 
the same problem. However, the exact embodiment of the multigrid paradigm depends 
strongly on the class of problems considered, with multigrid methods for differential 
equations (elliptic, parabolic, flow problems) being significantly different from meth- 
ods for integral equations. To place our problem in the context of multigrid methods 
we consider the simplified version of obtained by removing the inequality con- 
straints on the control, e.g., Uad = case in which reduces to the linear system 

(/C*/C + /3/)u = /C*yrf , (1.4) 

which represents the normal equations associated with the Tikhonov regularization 
of the ill-posed problem 

JCu = yd . (1.5) 

Beginning with the works of Hackbusch (see also [5]) much effort has been devoted 
to developing efficient multigrid methods for solving equations like (jl.4[) and (jl.Sp . 
e.g., see [ISifTBlfTIIlfHl fTllS] and the references therein. For example, Draganescu and 
Dupont [B] have constructed a multigrid preconditioner A4h which satisfies 

1-C— ^y— -^1 + C— , (1.6) 

13 {Mhu,u) (3 

where h is the mesh-size, Gh is the discretized version of {JC*JC + /3I), and p is the 
order of the discretization {p — 2 for piecewise linear finite elements). We regard (|1.6I) 
as optimal-order scalability since it implies that the condition number of the unpre- 
conditioned system cond(tJ;i) = 0(1//?) is been reduced by a factor of hP in the 
7W/i-preconditioned system, namely cond{Mj^^Gh) = 0[hP / (3); this reduction is of 
optimal order given that the discretization order is p. A similar situation is encoun- 
tered in classical multigrid for elliptic problems where multigrid is used to reduce 
the condition number from 0{h^^) to 0(1), the latter implying the desired mesh- 
independence property. We should point out that ()1.6|) implies that the number of 
iterations actually decreases with ft, J, to the point where, asymptotically, only one 
iteration is needed on very fine meshes. 

The presence of explicit box constraints on the controls and/or states in PDE- 
constrained optimization problems is sometimes critical both for practical (design 
constraints) and theoretical reasons (e.g., states representing densities or concentra- 
tions of substances have to be nonnegative). Methods for solving optimization prob- 
lems with inequality constraints are fundamentally different and more involved than 
those for unconstrained problems and generally fall into two competing categories: 
interior point methods (IPMs) and active-set methods such as semismooth Newton 
methods (SSNMs) . Both types of methods exhibit superlinear local convergence and 
can be formulated and analyzed both in finite dimensional spaces as well as in function 
spaces [171 [H E] i the latter being a critical step towards proving mesh- independence 
for the number of optimization steps. Both IPMs and SSNMs are iterative proce- 
dures that require the equivalent of a few PDE solves (i.e., applications of /C) for each 
iteration - called here outer iteration - and the solution of one or two inner linear 
systems at each outer iteration. Efficiency of the solution process is measured by the 
number of outer iterations needed to solve the problem to a desired tolerance (ideally 
mesh-independent) and by the ability to solve the inner linear systems efficiently. In 
this work we concentrate on the latter. 
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Even though they essentially solve the same problem, the linear algebra require- 
ments for IPMs are different from those of SSNMs. If formulated in the reduced 
form that is, with the PDE constraints eliminated (application of JC treated as 

a black-box) the structure of the systems arising in the IPM solution process is shown 
to be similar to the system (ll.4[) for the unconstrained problem i?,. More precisely, 
for IPMs we need to solve systems of the form 

gh,xu={K,*K, + Vx)u^b, (1.7) 

where I?a is the multiplication operator with a relatively smooth function A. Moreover, 
under specific conditions and for a natural discrete formulation of Draganescu 
and Petra have constructed multigrid preconditioners for the linear systems (|1.7p 
that exhibit a certain degree of optimality similar to the one multigrid preconditioners 
in [B]: the resulting multigrid preconditioner AAh for the Qt.x satisfies 

l-C^^I|A-^lk^^ j%^<l + C^||A-^|U^, (1.8) 
p {MhU, u) p °° 

We recognize in (jl.8|) the optimal-order h?-terva (linear splines were used for discretiza- 
tion), but also remark that the quality of the preconditioner normally is affected by 
the non-smoothness of A, which is generally expected to happen when the solution 
approaches the boundary. 

In this article we use similar ideas to design preconditioners for SSNMs. In this 
sense the present work should be regarded as a companion of [?]• We will show that 
the linear systems to be solved are essentially principal subsystems of (jl.4l) where the 
selected rows (and columns) correspond to the constraints that are deemed inactive 
at some point in the solution process. The constructed multigrid preconditioner is 
shown to essentially satisfy ()1.6|) with p = ^. While this order of approximation is 
clearly suboptimal, it still brings a significant reduction of the condition number if 
y/h <^ (3, and still results in a solution process that requires fewer and fewer iterations 
as h 10. 

We should point out that the strategies described above apply in general to the 
reduced problem (II. ip . However, a significant amount of literature is devoted to 
multigrid methods applied to the complementarity problem representing the KKT 
system of (|1.3p . Of these techniques we mention the collective smoothing strategy 
of Borzi and Kunisch [2l. For further references we refer the reader to the review 
article [3^ of Borzi and Schulz. 

This article has the following organization: in Section [2] we give the formal intro- 
duction of the problem, briefly discuss semismooth Newton methods, and present the 
main results. Section[3]is essentially devoted to proving the main result, Theorem l2.2l 
In Section 2] we show some numerical results to support our theoretical work, and we 
formulate some conclusions in Section [5l 

2. Problem formulation and main results. Our solution strategy will follow 
the discretize-then-optimize paradigm, where we first formulate a discrete optimiza- 
tion problem associated with (|l.ip . which we then solve using semismooth Newton 
methods. After introducing the discrete framework in Section lOl we give the details of 
the optimization method and its linear algebra requirements in Section[22I We define 
the two-grid preconditioner in Section 12.31 and state the main results. Furthermore, 
we discuss the multigrid preconditioner in Section [2.41 
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2.1. Notation and discrete problem formulation. Let i7 C M'' (d = 1, 2, 

or 3) be a bounded domain which, for simplicity, we assume to be polygonal (if 
d = 2) or polyhedral (for d = 3). We denote by W^ift), H"^{n), H^{n) (with 
p e [l,oo],m G N) the standard Sobolev spaces, and by || • || and (•, •) the L^-norm 
and inner product, respectively. Let H^™-{H) be the dual (with respect to the L^-inner 
product) of H"^{n) n i?o(ri) for m > 0, with the norm given by 

i^li?-'"(f2) = {u,v) /\\v\\Hm(_n^ . 

tie_f/™(o)n_ff,5(Ji) 

The space of bounded linear operators on a Banach space X is denoted by S-iX). 
We regard square n x n matrices as operators in £(R") and we write matrices and 
vectors using bold font. If A is a symmetric positive definite matrix, we denote by 
(u, v)^ = v'^Au the A-dot product of two vectors u, v, and by |u| a — \/ (u, u)^ the 
A-norm; if A = I we drop the subscript from the inner product and norm. The space 
of m X n matrices is denoted by Mmxn', if m — n we write M„ instead of Mmxn- 
Given some norm || • \\s on a vector space X, and T G S'iX), we denote by \\T\\s the 
induced operator-norm 

||T||, = sup \\Tul . 
uex, ||m|U=i 

Consequently, if T e £(L^(r2)) then ||T|| (no subscripts) is the operator-norm of T. 
If A" is a Hilbcrt space and T e £,{X) then T* e £,{X) denotes the adjoint of T. The 
defining elements of the discrete optimization problem are: the discrete analogues 
of /C, discrete norms, and discrete inequality constraints, all of which we introduce 
below. 

In the interest of the presentation we will make a few more specific choices for 
our optimization problem, namely we let W = 3^ = L'^{ fl), and reduce the constraints 
on the control to non-negativity, so that problem becomes 

minimize ^\\ICu — J/dlP + f /3 > fixed, 

subj. to : we L^{^), u ^ a.e. C^-^) 

To discretize the problem we consider a sequence of quasi- uniform (in the sense of [S]) 
meshes 7}, j = 0, 1, 2, . . . , which we assume to be either simplicial (triangular if d — 2, 
tetrahedral if d = 3) or rectangular, and let 

hj = max{diam(T) : T £ Tj} , j = 0, 1, 2, . . . . 

It is assumed that there are mesh-independent constants < / ^ / < 1 (usually 
1 = 7 = 1/2) so that 

/ < h,/hj^i ^ 7 . 
We define the standard finite element spaces: for simplicial elements let 

V^ = {ueC(n) -.yreTj, m|t is linear, u|ao = 0} , 
and for rectangular we use piecewise tensor-products of linear polynomials 



V^ = {ueC{n) : VTeTj, u\t e Qi, u\9n = 0} , 
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where 

f ' 

Qi < Cj Y\ lj.k{xk) '■ Ij.k linear polynomial of one variable 
[ J fc=i 

We assume that Tj+i is a refinement of Tj so the associated spaces are nested 

Since the algorithms and results are the same for both types of finite element spaces 
we will denote by Vj either V| or VJ. Let Nj — dim(Vj) and , . . . , P^J^ the nodes 
of Tj that lie in the interior of Jl, and define J} : C(0) — > Vj to be the standard 
interpolation operator 



(i) 

where f / , i — 1, . . . , Nj are the standard nodal basis functions. If we replace exact 
integration on an element T with vertices Pi, . . . ,Pi, by the cubature 

/ fi.)d. . ^ E /(P) , 
It 

'' ^ P vortex of T 

then the L^-inner product is approximated by the mesh-dependent inner product 
{u, v)j = J2 ^P'' u{pP)v{P^'^), for u,v€ Vj , 



where 



E (2.2) 

The discrete norms are then given by 



lu\\j =^ sj {u, u). . 

Since the quadrature/cubature is exact for linear functions, or tensor-products of 
linear functions, respectively, we have 

{u,v)- — / Jj{uv), for all u,w S Vj . 



n 

Moreover, due to quasi- uniformity, there exist positive constants Ci,C2 independent 
of ^ such that 

CiM |||«|||, s;C2||?/||, Vug V, . (2.3) 

We should point out that the norm-equivalence (|2.3p extends to show mesh-independent 
equivalence of the associated operator-norms. We say that the weights wp'* are uni- 
form with respect to the mesh Tj if there exists Wj > independent of i so that 

u)p^ ^oj^h'^ for i^l,...,Ni . 
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We call a triangulation locally symmetric if for every vertex P the associated nodal 
basis function ip is symmetric with respect to the reflection in P, that is, 

ip{2P -x) = ip{x), yx en . 

(i) 

If a mesh is uniform then it is locally symmetric and the weights are uniform. 

On each space Vj consider an operator ICj G S^iVj) representing a discretization 
of /C. For the discrete operators we denote /C* to be the adjoint of ICj with respect to 
(•, ■)-, that is, 

{lC*u,v)_. = {u,ICjv)^ , yu,v e Vj . 

We assume that the operators satisfy the following condition: 

Condition 2.1. There exists a constant C — C{1C) depending on IC,fl,To and 
independent of j so that the following hold: 

[a] smoothing: 

||/Cu|U,„(a) C \\u\\ , Vu e L2(r!), m = 0, 1, 2 ; (2.4) 

[b] smoothed approximation: 

WKu - ICjulH^^n) ^ Ch]--^ \\u\\ , -iu e Vj, m - 0, 1, j ^ ; (2.5) 

[c] uniform boundedness of discrete operators and their adjoints: 

max(||/C*M 1^00(0), ||/C-,-u||L~(n)) < C ||u|| , Vu e V^, j > . (2.6) 

We now formulate the discrete optimization problem using the discrete norms 
and non-negativity at the vertices: 

minimize i|||/CjU - y4j + §Mj, /3 > fixed, 

t-^ (2-7) 
subj. to : ueVj, u{P>^') ^ for i = 1, . . . , iVj . 

The formulation (|2.7|) is identical to the one described in [7], except for we retained 
only one inequality constraint for simplicity. It is worth noting that the pointwise 
imposed non-negativity constraints at the nodes imply the non-negativity of the con- 
trol everywhere. This is not the case with higher order elements where nodal basis 
functions can have negative values. 

If Kj is the matrix representation of K,j in the nodal basis and Wj is the diagonal 
matrix with diagonal entries w^"''' , ^2"'^ , • ■ • , w'fj^ , then (j2.7|) reads in matrix form 

minimize J/3(u) =^ i |KjU - yd|^. + | |u|^. , /3 > fixed, 

' " ' (2.8) 

subj. to : u e M^J , u ^ . 

Furthermore, we write Jfj{u) — ^u^CjU — bju + 7_,-, where 

C,=KjW,K,+/?W„ b,=KjW,yd, 7, = ^ylW.y, . (2.9) 

We also point out that the adjoint operator /C* is represented by the matrix Wj^KjWjKj 
Throughout the remainder of this article we will omit the subscript j when focusing 
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on a single grid. Since Jp is strictly convex and quadratic, (j2.8l) has a unique solution 
u, which is given by the Karush-Kuhn- Tucker (KKT) conditions: there exists a vector 
A e M.^^ so that 



u A = , 

u ^ 0, A ^ , 



A = , 



(2.10) 



where u • A is the vector [uiAi, U2A2, ...]-'". 

2.2. Semismooth Newton methods. To ease notation, throughout the first 
part of this subsection we leave out the subscripts "j" if there is no chance of confusion. 
The KKT complementarity problem (I2.10p can be formulated as 



Cu b - A = , 

A — max(0, A — /3u) 







(2.11) 



As shown by Hintermiiller and Ulbrich in [12] , the semismooth nonlinear system p. lip 
can be solved very efficiently using semismooth Newton methods, also shown in |11] 
to be equivalent to the primal-dual active set method which we now describe briefly. 
An immediate consequence of the second equation in (|2.11[) is that A, u ^ 0. A second 
consequence is that if for some index i we have A^ — j3\ii > then = 0, so the 
i"^ constraint is active. Instead, if A.; — /3ui ^ then A^ — 0, and the constraint is 
deemed inactive. Hence we define the active index-set by 

A = {ie{l,...,N,} : (A-/3u), >0} 

and the inactive index-set by 

I^{ie{l,...,N,} : (A-/3u), «;0} . 

The semismooth Newton method produces a sequence of active/inactive sets {Ak,Ik)k=i.2,. 
that approximate {A,T). Given (AkT^k), we set the system 



r cuC^+i) - a(*^+i) = b , 

('=^+1) = 0, for i e Ak 
0, for i elk ■ 



(2.12) 



Ak+l) 



We divide the vectors u ~ u^'^^^^ and A = A^'^^^-' into their active and inactive com- 
ponents, and we also block-divide the matrix C accordingly: let C^^ — {Crs)r,seXk be 
the principal minor associated with the inactive constraints, C^'^ = {Crs)reXk,seAk: 
etc. With this notation, (|2.12p reads 








(jAA 
I 





I 




I 





I 






the solution of which is 









" hi ' 










hA 






A/ 







1 













, 


and 




-- C^'ui - 



(2.13) 
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The new active and inactive sets are given by 

Ak+i={ie{l,...,N,} : (A('=+i)-/3u(^+i)), >0} 
Ifc+i = {*€{!,..., iV,} : (A('=+i)-/3u('=+i)), SCO} . 

Hence the critical system to be solved is at each semismooth Newton iterate is 

= . (2.14) 

After dividing the i*^ equation in (|2.14p by w^'''' we obtain the system 

Gfuf+'^^{Wj%)i. (2.15) 

where G^" is the principal minor of the matrix (W^^KjWjKj + /3I) corresponding 
to the inactive constraints indexed by Ik- We should point out that for large-scale 
problems the matrices Kj are expected to be dense and are not formed. Thus we 
can solve (|2.15p only by means of iterative solvers, and therefore efficient, matrix-free 
preconditioners are necessary for accelerating the solution process. Also note that if 
Ak = 0, then the multigrid preconditioners developed in [5] can be used. Our main 
contribution in this work is the design of two- and multigrid preconditioners for (j2.15p 
for the case when Ak 7^ 0- 

A final remark to related to the system (I2.15|) refers to the use of the diagonal 
matrices Wj. These should be viewed as diagonal mass matrices obtained via the 
discrete inner products (•, If we used the exact L^-norms instead of the mesh- 
dependent norms in (12.71) . then (|2.14p would still be the correct system to solve, but 
it would not be reduced to (|2.15p . As will be shown in the Section 12. 3[ it is precisely 
the form (|2.15l) that lends itself to efficient multigrid preconditioning. 

2.3. The tviro-grid preconditioner. Let j ^ 1 be a fixed level, to which we 
shall refer as the fine level. We will first define a two-grid preconditioner for that 
will involve inverting Gj'Li- To design the two-grid preconditioner we will regard the 
matrix G^" as an operator between finite element spaces. 

2.3.1. Construction of two-grid preconditioner. Let Z^^' C {1, . . . ,Nj} be 
the set of indices corresponding to inactive constraints at some outer iterate considered 
fixed (this set will change at each outer iteration). We will call a node or a vertex 
inactive if the corresponding constraint is inactive. Define the fine inactive space by 

Vf =span{^p) : lel'^^^} , 
and the fine inactive domain 

y supp(^p)) , (2.16) 

where supp(u) is the support of the function u. The critical component of the pre- 
conditioner is the definition of the coarse inactive index set: 

lO-i) {l,...,A^,-i} : suMvf') C . (2.17) 

If if is the index in the fine numbering associated to the coarse index ic, that is, 
pj;p = pj;^ ^\ the definition above is equivalent to saying that ic G l'-'^^^ if i/ 
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together with all its neighbouring fine indices are inactive. In Figure \TJ\ we depict a 
set of inactive fine nodes by filled circles, and the associated coarse inactive nodes by 
hollow circles. The coarse nodes that are not inactive are shown with a square hollow 
marker. The coarse inactive space is now defined to be 

Vj-ii = span{(^p"') : i G I'^J-^^ , (2.18) 

and the coarse inactive domain is now 

nf^,= U supp(v.?-^V 

Note that Vj" C Vj" and C f2™. In connection with ft"^ we also define numerical 
interior Intufi^" of f2^" (relative to ^]jLi) to be the union of all coarse elements T 
included in the fine inactive set, that is, T C $7-^", and whose vertices are either in 
or lie on the boundary of fl (see Figure [^TT|) . Furthermore, let the numerical 
boundary of il™ (relative to ^"-i) be given by 

Note that Intn^j" C flf_^. 

Let now Wjl^ be the L^-orthogonal complement of Vj^i in Vj" and define the 
L^-projectors 

^in . -\ liii , -\;in ^in . -\ liii , ■\y\;in 

SO T^Y-i + pT-^ identity on V™. Furthermore, let £™ : V™ —?> Vj be the natural 

embedding obtained by extending a function with zero outside of We may oc- 
casionally skip the explicit use of the embedding operators to ease notation. We also 
define the restriction Rj-i : Vj — > Vj-i as the adjoint with respect to to the 

embedding of Vj-i in Vj, that is, for u E Vj 

{u, v)^ - {Rj-iu, v)^_^ , Vf G V,_i , (2.19) 
and let Pj" : Vj — > Vj" the projection with respect to (•, ■) ■ given by 

\i=l / iGlO) 

Furthermore, denote by Hj the orthogonal L^-projcctor onto the space Vj. From the 
equivalence (|2.3p of the discrete norms with the i^-norm it follows that 

\\R,u\\^C\\u\\ , J = 0,1,..., . (2.20) 

for some mesh-independent constant C. 

The matrix G^" in (|2.15l) represents the operator 

gf - Pf (/C*/C, + /?/) £f , (2.21) 
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Fig. 2.1. Inactive fine nodes are marked with a dot, inactive coarse nodes are marked with a 
circle, and active coarse nodes are marked with a square. Note that the coarse nodes A, B, C are 
active because the supports of the corresponding nodal basis functions are not included in f2^° , even 
though the nodes themselves lie in the interior of OJ" . The area between the solid and the dotted 
lines lies in the numerical boundary of H'." . 



and thus matrix-vector products for G" are computed accordingly. We define a 
two-grid preconditioner il/j" for Q™ as in Draganescu and Dupont for the 
unconstrained case: 



Mf - PfL^ (/C*_i/C,_i + pi) ^j-Li + . (2.22) 
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Note that the inverse 5" of the preconditioner Afj" has the exphcit form 

^„ drf (^j^inyl ^ (gin ^fl^ + pf_^ . (2.23) 

Since none of the matrices representing C/j" or A^" are formed, it is the operator Sj^ 
that we need to apply in practice, so the explicit formula p.23p for (A^^")~^ is es- 
sential. In light of this fact we should remark that, if the projection Trjl^ were to 
be replaced by a restriction operator, as is the case in classical multigrid, then ()2.23p 
would no longer be true. 

2.3.2. Matrix-form of the preconditioner. In order to describe the matrix- 
form of the preconditioner in Matlab form we introduce the mass matrix Lj on Vj , 
and we denote by Jj S Mn-xn^^i the matrix representing the interpolation operator 
Jj\vj-i G £(Vj_i, Vj). Furthermore, let Rj-i = 2~'^Jj S Mn^_ixNj be the restric- 
tion operator. We assume that the inactive indices on level j are stored in the vector 
ij, and define the matrices 

L}=L,(i„i,), Jj=j(i„i,_i), R}_i=R(i,_i,i,), E} = I(:,i,), P} = EJ, 

where we used Matlab notation for the selection of submatrices. Note that Ej rep- 
resents the extension operator £™ and Pj the operator P™. We can now write the 
projector-operator t^^-i matrix- form as 

and p^j^-i is represented by (I — Jjnj_]^). So A^^" is represented by the matrix 

represents ^ 

Mj = J} Pj_i(Kj_iK -f /3I)Ej_i nj_i + p{i - Jjnj.i) , 
and 5j" is represented by 

s} = J} (p}_i(Kj_iK,_i + /3i)Ej_0"' nj_i + r'ii - Jln}_i) . 

We should point out that, due to the presence of n^_i, the matrices Mj and Sj are 
slightly nonsymmetric, hence one has to employ solvers for nonsymmetric systems in 
connection with the preconditioner. We found that conjugate gradient squared 
(CGS) is quite efficient (see Section 

2.3.3. Spectral distance estimation. To quantify the quality of the two-grid 
preconditioner A^^" we will estimate the spectral distance da between {Q™^)~^ and 
its inverse 5™ given by (|2.23p . The use of the spectral distance ensures that such 
estimates extend automatically to multigrid preconditioners, as shown in [51 [7] (see 
Section We briefly recall the definition of the spectral distance, as introduced 
in [B]. Given a Hilbert space (A", (•, •)) we denote by 2,+ {X) the set of operators with 
positive definite symmetric part: 

Z+{X) = {T e 2,{X) : {Tu, u) > 0, e X \ {0}} . 

Let the joined numerical range of S*, T e £+(<%") be given by 

I (Tcw, w) 
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where Tc{u + iv) — T{u) + iT{v) is the complexification of T. The spectral distance 
between S,T Cz is a measure of spectral equivalence between S and T, and it 

is defined by 

d^iS,T) = sup{|lnz| : z € W{S,T)} , 

where In is the branch of the logarithm corresponding to C \ (— cx),0]. Following 
Lemma 3.2 in fS', if 1^(5", T) C Sa(l) = {z e C : |z - 1| < a} with a e (0, 1), then 

d„(5,T) Mii:Msup{|z-l| : zeW{S,T)}, (2.24) 
a 

which offers a practical way to estimate the spectral distance when it is small. The 
spectral distance serves both as a means to quantify the quality of a preconditioner 
and also as a convenient analysis tool for multigrid algorithms. Essentially, if two 
operators S, T satisfy 



1-6^ 



{Scw,w) 



{Tew, 1 



with (5^1, then dcr{S,T) Ri 6. If « G ^ is a preconditioner for G, then both 
da{N,G-^) and da{N-^,G) (quantities which are are equal if G, N are symmetric) 
are shown to control the spectral radius p{I — NG) (see [3) which is an accepted 
quality-measure for a preconditioner. The advantage of using d^ over p{I — NG) is 
that the former is a true distance function. 
The main result of this article is 

Theorem 2.2. // the operators K, and ICj satisfy Condition 12.11 and the weights 
are uniform, then there exists S > and a constant G(/C) (see Condition 12.11 ) 
independent of j and the inactive set so that 

d, ({Qfy\sf) «c G/3-1 [h] + , (2.25) 

where fi™ is the Lebesgue measure of d-oTt™ , provided that 

(/i- + v^) < s . 

We postpone the proof of Theorem 12.21 until Section |3l 

Remark 2.3. The natural question arises as to how to estimate /i™. In the worst 
case scenario there are no coarse inactive nodes, so = 0; therefore dn^lj^ = ^, 

case in which the two-grid preconditioner is 131 , so essentially there is no precondi- 
tioner. However, if u is the solution of (|2.ip and the continuous inactive set defined 
by ri'" = {a; G : u{x) > 0} is a domain with Lipschitz boundary, then the discrete 
inactive set ^l™ is expected to be close to 57'" provided that a good initial guess at the 
inactive set is available. In this case we expect that d^Vt™ will lie within Ghj of the 
topological boundary offl™, therefore 



fif « Ch 



1 ' 



where G is proportional to the {d — 1) -dimensional measure of . Hence the esti- 
mate (|2.25p truly implies 



{{gf 



)-\sf)^G^ , (2.26) 
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which is consistent with the numerical experiments in Section^ 

In case the grid is quasi-uniform but not uniform we apply Theorem 12.21 to the 

— dcf 

matrix Kj = WjKj. The important aspect in the estimate is the verification of 
Condition 12.11 by Kj . Fohowing [7] we introduce the following indirect measure of 
grid-smoothness: for each j we consider a C^-function : O — >■ R so that 

w,{pI'^)^wI^\ Vz = l,2,...,iV,, (2.27) 

with wp-* given by (|2.2p . With this notation the matrix Kj represents the opera- 
tor JCj e £(Vj) defined by 

ICjU ^'^ Jj {wj ■ (JCju)) . (2.28) 

Note that, because the grids are hierarchical (7j is obtained from Tj-i by adding 
nodes), the function Wj can serve for defining all operators ICi, for I = 0,1,..., j. 
Consider j fixed, and define JC E £(L^(f])) by 

dcf / 1^ \ 

Itu = Wj ■ (K.jU) . 

By Proposition 4.8 in [7j, the operators {ICi)i=o^i^,,,j together with the continuous 
operator /C satisfy Condition 12.11 with 

C(^) = lk. lkr(n)C(/C) . 

Thus we establish the following 

Corollary 2.4. If the operators JC andJCj satisfy Condition 12.11 and Wj e C^(f7) 
satisfies ^2.21^ , then there exists (5 > and a constant C(1C) independent of j and the 
inactive set so that 

da ({Qfy',Sf) ^ Cr'WwjWw^^i^) (h, + , (2.29) 

where /i" is as in Theorem \2.2[ provided that 

We should remark that the power of hj is 1 as opposed to 2 in Theorem 12.21 This 
fact is due to the nonuniformity of the grid as can be seen in the analysis of Section [3] 

However, in general the larger of the two terms in (|2.29p is y^Mj^, so the main effect 

of the nonuniformity on the estimate is the presence of the factor ||'u;j||vK2°°(o)- 

2.4. The multigrid preconditioner. We now assume the levels j — 1 ^ jo ^ 
to he fixed (we refer to jo as the base- level), and the goal is to construct a multigrid 
operator so that -Ej"^ « (^j") which satisfies the following conditions: (i) 

Zf^^^ = Sf; (ii) the estimate (p?^ holds if we replace Sf" with Zj^^^. In order to 
construct -Zj"^ we must first specify, starting at the finest level j, the coarser inactive 
domains fi™, inactive index-sets l'-*^-' , and inactive spaces V^" for fc = j — 1, . . . , jo. All 
these entities are defined recursively using (|2.16p . ()2.17p . and ()2.18p and are essentially 
specified by the sets Z''^^. Hence we give below the algorithm for computing the 
inactive-index sets I*^'^^ for fc = j — 1, . . . , jo (note that X^^^ is given by the semismooth 
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Newton method iteration). Given a vertex pj^'^ of the triangulation Tk-i, let P/*^-* 

(k—l) 

be its fine labeL we define the "fine neighbourhood" of P^^ by 

J^kiP-^'^^) = {R^'''> ■■ R^'''' neighbour in Tk of pj;^^} U {pj;^^} . 

Algorithm 2.5 ( (Inactive set definition)). 

1. for k = j : -1 : jo + l 

2. xC^-i) = 

3. f or i = 1 : Nk-i 

I if A4(/^^'"'^) ClW 

5. j{k-i) ^J{k-l)^J^p{k-l)y 

We denote by tt^" : L^(J7) V^" the i^-projection and we define the operator 

3ti : ii(Vti) ^ £(Vr), atil-^) = A- . + /3-i(/ - ^iP) . 

We should point out that the operator can be written as 

Sf^ = {{Gr-,)-') . (2.30) 

In fight of equahty (I2.30p and the continuity of the affine operator it is tempting 
to define the foUowing muhigrid preconditioner: 

(gf)-l , if J = JO, 

. . (2-31) 

3i-i(2j-ijo) , ifj-l^jo- 

As pointed out in the T^-cycle type preconditioner -Zj"^ does not satisfy condi- 
tion (ii) above. In fact one can see numericaUy that, under the conditions set in 
Remark 12.31 satisfies (|2.26p with hj replaced by Hq. As a result the number of 

preconditioned iterations would no longer be decreasing with hj \. 0, as is the case for 
the two-grid preconditioner, but would be fixed. Instead of the definition (j2.3ip . we 
employ the same strategy adopted in |6l [7] , which guarantees that the estimate for 
the multigrid preconditioner will essentially be the same as the one for the two-grid 
preconditioner (except for a constant factor) . In order to do so we define the operator 
9^fc by 

Dlfc : £(V[") ^ £(Vi"), mk{X) '^^^ 2X-X -Gi^'-X . 
Note that ^k is the Newton iterator for the operator-equation 

- GT = , 

as shown by Draganescu and Dupont in [B]. This implies that, if X^ is a good approxi- 
mation of {G]^)~^, then Xi = OTfc(A'o) is the first Newton iterate of the above operator- 
equation starting with Xq, and so Xi is significantly closer to (CJ™)~^ than Xq. This 
idea was also used in '7 to construct multigrid preconditioners of the same quality as 
the two-grid preconditioners. 
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Algorithm 2.6 ( (Operator-form definition of Zj"^-^)). 

1. if in =i-i 

2. Zfj^ = Sf 7. two grids 

3. else (here we expect j — 2 ^ jo) 

I Zf:j„ = I multiple grids 

Since the behaviour the two-grid preconditioner is less than the optimal, as noted 
in Remark I2.3[ we choose not to present an ellaborate analysis for the multigrid 
preconditioner; suffice it to say that, under the setup of Remark 12.31 and by using 
similar heuristical arguments together with the formal arguments in [B], we argue 
that the following estimate holds 

provided the coarsest mesh-size /ij^ is sufficiently fine. Naturally, the cost of applying 
Zj"^ is lowest if jo is minimized, yet a larger jo ensures that Zj" ^ has an efficiency 
that is comparable to S™. However, for truly large-scale applications, e.g., for c? = 3 
or 4, a choice of jo = j ~ 2 brings a significant reduction in size for the coarsest 
problem and may be sufficient. 

3. Analysis. The main step in the analysis is to evaluate the norm-distance 
between the operators Q™ and AA™ which is done in Proposition 13.41 The plan of 
the analysis generally resembles that of the analysis of multigrid preconditioners for 
interior point methods from [i7|, however, certain critical estimates related to the 
projection 7rj" are different for the case of semismooth Newton methods. 

First we restate Lemma 4.3 in [7] as 

Lemma 3.1. With{w'f ^) i^i^N chosen as in p.2p there exists a constant C ~ C(7o) > 
independent of j so that 

\{u,v)^ - {u,v)\ < Ch'^j\\u\\Hi(n) ■ Mm(n) , ^u,v £ Vj. (3.1) 

We also recall Lemma 4.4 in |7 : 

Lemma 3.2. If ]C,ICj satisfy Condition \2.1\ there exist constants C{K.) and 
C = C"(r2) independent of j such that the following hold: 

(a) H^,L^ - uniform stability oflCj: 

l|/C,"|U..(o) ^ CilC)\\u\\ , e V„ m = 0, 1, j = 0, 1, . . . ; (3.2) 

(b) smoothing of negative-index norm: 

\\ICu\\ s$ C(/C) ||u||jj-,„ , V-u e Vj,TO = 1,2 ; (3.3) 

(c) negative-index norm approximation of the identity by Hj^i^TZj^i: 

1(7 - ^,-i)w||5-.(^) ^ C'/il , e V,; (3.4) 
\\{I-n,^i)uy_,^^^^C'hP\\u\\, VueV,, (3.5) 

where p — I on a quasi-uniform grid, and p — 2 on a locally symmetric grid; 

(d) K. diminishes high-frequencies: 

||/C(/ - TTj-i)u\\ ^ C{K)h] \\u\\ , Vu e V, ; (3.6) 
\\JC{I - 7^J-l)u|| C(/C)/i^ ||u|| , Vu e Vj , (3.7) 
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where p — 1 on an unstructured grid, and p = 2 on a locally symmetric grid; 
(e) 

\{lCu,JCv) — {JCjU,JCjv)\ ^ C(/C)/i|||u|| • ||ti|| , Vu,w e Vj . 



(3.8) 



The main difference between the two-grid preconditioners for the unconstrained case 
versus the constrained case is in the properties of the projectors on the coarse spaces. 
While the projector ttj-i on the entire coarse space Vj-i satisfies p.4p . the projector 
on the inactive space Vj" satisfies the weaker estimate below. 

Lemma 3.3. If d ^ 3 there exists a constant C depending on the domain ^l and 
the base triangulation To so that 



\\ul for aUwe vj" , 



(3.9) 



where ij}^^ is the Lebesgue measure of 8^^™. The constant C is independent of j and 
of the inactive set flj^. 

Proof Let u G Vj", v e H'^{n) f] H^{n) be arbitrary, and let Vj-i G Vjl^ be the 
natural interpolant of v in Vj"_]^, that is 



Let T e Tj-i be a coarse element lying in ft"^. If T C Intn^^j" then Vj-i agrees 
on T with the interpolant of v in Vj_i. Therefore a standard interpolation estimate 
(see [S]) applied on Int„r2^" gives 

Ik- Vj-i|lL2(int„a'») Ch^j^^\v\m{int„np ^ Cr^h^j\v\H2(n) . 

On a coarse element T that satisfies T C dn^]p, if Vj-i is not identically zero on T then 
Vj-i and V agree at least for one vertex of T and potentially disagree at a vertex cor- 
responding to an active constraint. In either case the bound ||l°°(t) ^ lkllL°°(T) 
holds. Hence 



by Sobolev's inequality. We have 



= \{{I-7tJL,)u,v-v,^,)\ = 



Int,,f2" 



(u - 7r'"_iu) {v - Vj-i) 



< C\\u - 4liu|L.(fji„) + y;^) \\v\\h2 



(3.10) 
(3.11) 

(3.12) 
(3.13) 



Since ||7rj"_2''^ll ^ I""!' the result now follows after dividing by ||w||//2(f2) and taking the 
supremum over all v e H^{n). □ 
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Proposition 3.4. // the operators K, and ICj satisfy Condition 12.11 and the 

weights wf are uniform, then there exists a constant C independent on j and the 
inactive set so that 



\\gr~Mf\\^c(^h'^ 



(3.14) 



where /i™ is the Lebesgue measure of d^i^J^. 

Proof Since = nfl^ + pf_-^ = /yn, we have 



(3.15) 



The argumentation has a structure similar to the proof of Proposition 4.5 in [7] with 
changes due to the specific approximation properties of PfLi and '^"-i- Let u,v € V™ 
be arbitrary, and define u = £jLiT^"-iU. Note that ||u|| ^ ||u||. We first examine the 
difference 



{I ~ P^l,)ilCUIC,-iu) ■ V 



(I - P]l,)ilCUIC,-iu) ■ V 



d„n" 



^ ||(/-P;'ii)(/C*_i/C,_iSf)|L.(9„o.") • |klL.(a„o.") 



where we used that {I — Pj"_i){IC*_iJCj-iu) is zero in Int„f2™ and that v is zero outside 
of O^". Since for w £ and x E ft the function value (/ — Pj° ]^)w(x) lies between 
and w{x) we have 



!(/- p-i)(/c;_i/c,_iS)|L.(a„o.=) ||/c*_i/c,_i^i|U^(a,) < c||s;|| < c\\i 



Therefore 



\{lC*_,IC,^,u,v)^ ~ {P^,iq_,IC,-iu,v)\ ^ cJtif \\u\\ ■ \\v\\ 



(3.16) 
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We return to estimating {Q"^ — A^^"). We have 



1239 



At 



Ao + {Kj£fu,K.3£fv) . - {K.j£fu,Kj£fv) + 



An 



{lCj£fu, K.j£fv) - {K.£fu, K.£fv) + 

A. 



{lC£fu,K.£fv) - {lC£f_iTTf_iU,IC£]\Rj-iv) 



{lC£ft,7Tfl,u, IC£fURj-iv) - {lCj-i£fUTrf_^u, K.j-i£f_^R,-iv) 

As 



(jC,-i£^^,u,K.,-i£JLiR,-iv) ~ {lC,-i£fL,7TjL,u,IC,-i£fL,R,-iv)^_^ . 
By (1311)) 

\Ao\^C^\\ul\\v\\ . 
For Ai and A5 we use (jSH]), (|3^ . and (|2?20| to conclude that 

l^il s; C/i2||/c,£f 7.1^1(0) • ||/C,ff ^ . ||„|| , 

\A,\ ^ Ch]_,\\JC,^,£fL,nf_,u\\H^n) ■ \\JC,£f^R,-iv\\H^n) 
^ Cfh',h^Uu\\ ■ \\R,-iv\\ < Ch',\\u\\ ■ \\v\\ . 

Estimation of ^2,^4, essentially involving p. 81) . is handled exactly as in 7] to give 

mSix{\A2\,\A4\)^Ch^j\\ul\\v\\ . 
Finally the term A3 is estimated by 

\A,\^\{JC£^{I -nJL,)u,JC£^v)\ + \{JC£fL,7r];^,uAIC£^{I ~ R,-i)^^ 

< c (/i^ + .y;;;^ ) 1^,1 . i^;! . 

We conclude that 

{{gf - Mf)u,v)^\ ^ C [h] + y^) \\u\\ ■ \\vl yu^ve , 



hence |||^™ — A^j"|||j ^ (jij + \J y^f^ , and the result follows from the equivalence of 
the operator-norms ||| • |||j, || • ||. □ Theorem l2.2l now follows easily from Proposition 
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Proof. First note that for u e Vj" 

(gf u, u)^ - {lCf£"^u, ICf£"^u)^ + p {u, u)^ ^ p {u, u)^ , (3.17) 

so (j{Q"^) C [(3,oo), which imphes, due to the symmetry of Q™ with respect to (•, ■) ^ 
and the equivalence of the operator-norms ||| • |||j, || • ||, that 

II {Gfy^-i^cr'- ■ 

By Proposition 13.41 

||/ - [gfy'^ Mf (gf)'^ II ^ c|| (g]")"^ p • i^f - x;"|| ci3-^ [h] + . 

The remainder of the argument proceeds as in the proof of Theorem 4.9 in [J. □ 

4. Numerical experiments. We test our algorithm on the standard linear 
elliptic-constrained distributed optimal control problem 

minimize \\y-ydt + ^\ut ^^^^ 

subj. to —Ay = u , y £ H^ijl), u ^ a.e. in 57. 

where A is the Laplace operator acting on i?Q (il), so /C = (— A)~^. The problem (|4.ip 
is often used as a test example for multigrid algorithms (e.g., see [SI [7]). Parts [a, 
b] of Condition 12.11 follow from standard estimates for finite element solutions of 
elliptic equations, while condition [c] follows from standard L°°-estimates [4]. Note 
that the L°°-stability in [c] does not require optimal i°°-convergence rate of the 
discrete solutions to their continuous counterparts, but merely convergence 

lim \\ICjU — /Ciilioo — . 

We test our multigrid preconditioner in two contexts. First we consider an "in vitro" 
experiment where we artificially fix a nontrivial subset of the domain n based on 
which we construct (artificial) inactive sets on all grids. The goal of this experiment 
is to assess the change in quality of the two-grid preconditioner while varying only 
the mesh-size. For each grid we then construct the two-grid preconditioner and we 
estimate numerically the decay rate of the spectral distance da{{Gj^)^^ ,S"^). The 
advantage of this approach is that we can isolate the effect of the inactive set on 
the quality of the preconditioner, since in the context of actually applying the semis- 
mooth Newton method the inactive sets are expected to change from one resolution 
to another. The second test is an actual, "in vivo" application of the multigrid pre- 
conditioner in the context of solving (j4.1l) using the semismooth Newton method. 

4.1. One-dimensional "in vitro" experiments. Let /3 = 1, f2 = [0, 1] and 

designate the interval il™ = [1/8,3/4] as the set where "inactive vertices" reside; this 
set is used for all grids. Let rij — 16 • 2^ ,j ~ 0, 1, 2, ... , and define Tj to be the 
uniform grid with Uj intervals on 51 with mesh-size hj — 1/ nj . The "inactive indices" 
are given by 

2]" = {^ e {l,...,nj- - 1 : i h, e ff"} . 
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Table 4.1 

Spectral distance decay for a fixed inactive domain (fS = I). 



grid-size (rij) 


16 


32 


64 


128 


256 


512 




0.0023 


0.0016 


0.0011 


7.4617- lO-'' 


5.1996-10-'' 


3.6372 - 10-4 


dj-i/dj 




1.4610 


1.4506 


1.4421 


1.4351 


1.4295 



Correspondingly, we construct the matrices G^" and M^" as in Section 12.3.21 and we 
compute the quantities 

dj- = max{|lnA| : A G a(Gf , Mf )} , j 1, 2, . . . , 

where cr(A, B) is the set of generahzed eigenvalues of the matrices A, B. In general 
we have 

but if both matrices are symmetric then the above inequality becomes an equality. In 
this case G™ is symmetric and M™ is close to being symmetric, so we expect that dj 
is a good approximation of ^^(G"^ MJ"). We report in Table HTT] both the numbers 
dj and their ratios dj-i/dj for j = 2, ... ,6. The last row of the table indicates that 
the numerical results are consistent with the estimate dj ^ C ^Jh] in Remark 12.31 
namely we notice that 

dj-i/dj V2 . 

Another advantage of the one-dimensional example is that we can easily compute 
numerically all the generalized eigenvalues of G" and M^". What may be surprising, 
is that, while the largest of them decay a the predicted rate of ^yhj as hj | 0, for 
each grid most of the generalized eigenvalues are very close to 1. In fact, for a fixed 
grid, if we order the generalized eigenvalues according to their distance from 1, we 
notice an exponential decay of |Ai — 1| (see Figure |4T|) . We also plot in Figure |4?2] the 
function u G V2" associated with the generalized eigenvector u that corresponds to 
the generalized eigenvalue which is furthest from 1. What we notice is that u is large 
around the boundary of f2™. These facts suggest that, with the exception of a limited 
number of generalized eigenvalues associated with eigenvectors strongly related to the 
boundary of fi'", the generalized eigenvalues G™ and M™ are in fact quite close 1. 
We revisit this idea in Section 14.31 

4.2. Tvifo-dimensional "in-vivo" experiments.. We now consider the two- 
dimensional version of ()4.1|) with = [0,1]^. To construct uniform triangular grids 
7}, j = 0, 1, ... , we first divide each side of uniformly in Uj —64-2^ intervals to 
obtain a uniform rectangular grid, and we further divide each resulting grid-square 
in two triangles along the diagonal of slope —1; the mesh-size is thus hj = l/fij, and 
the number of variables is Nj = [rij — 1)^. The Poisson problem is then discretized 
using standard continuous piecewise linear elements on a triangular grid. We solve 
the problem using grid-sequencing, that is, the solution on level (j — 1) is used as an 
initial guess for the semismooth Newton iteration at level j. More precisely, if 
denotes the coarse inactive domain as defined in Section [273l we define the initial guess 
at the inactive set on level j by 

2^-')={ze{l,...,iV,} : supp(^«) C f];- J . 
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, a--B--B--i>-H--H--a--Q- - a- -B - - Q- - a- -□ 



Fig. 4.1. Absolute values of the generalized eigenvalues o/ G!," and Mj" (n2 = 32). There are 
21 eigenvalues, and they decay rapidly to 1. 




Fig. 4.2. The generalized eigenvector corresponding to the generalized eigenvalue Glj" and M!," 
which is furthest from 1. 



Then we apply the semismooth Newton iteration described in Section 12.21 and we 
solve the reduced linear systems (j2.14p at each outer iteration using the conjugate 
gradient squared (CGS) method preconditioned by the two-grid explicit precondi- 
tioner Mj^- For comparison we also solve the same systems using unpreconditioned 
conjugate gradient (CG). Since the process is matrix- free (we solve the Poisson prob- 
lem using classical multigrid) we use the explicit form S"^ of the inverse of A4"^ as 
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given in (j2.23p . For each outer iteration we record the number of preconditioned CGS 
iterations required for solving the hnear system (I2.14p . For the numerical example we 
choose as "target" control the function defined by 



Udix) 




if \\x - xo\\ < r 
otherwise 



for some xq E fl (see Figure 14. 3p , and let y be the solution of the Poisson equation 

-Ay ^ Ud, v\dn = . 

The "data" yd entering the control problem (|4.ip is obtained by adding a random 
perturbation to y, namely 

Vd^y + S, where \5\l-^ ^ 0.05 ||y||L=° • 

If 0, and if 5 and (3 are small then the solution w™'" of (|4.ip is expected to be 
close to Ud- Therefore, by letting the parameter a to be slightly negative we expect to 
find a localized u™'" with nontrivial active and inactive sets, as desired for testing the 
performance of our algorithm. In Figure we show the solution u™™ for a = —0.1 
and TiQ — 64; as it turns out, for this example about 11.5% of constraints are inactive 
if /? = 10^^, and about 51% are inactive if /? = 10^"'. 




Fig. 4.3. Target control with zq = [0.54,0.62]^, r = 0.06, a = -0.1. 

The data in Table l4?2l shows the iteration count for the two-level preconditioned 
CGS method for (3 — 10^'* and for unprcconditioned CG. For example, on the third 
grid (n2 = 256) the semismooth Newton method converged in 4 iterations, the first 
converging after 4 preconditioned CGS iterations T(or 13 CG iterations), while the 
remaining three outer iterations required 5 CGS iterations each (or 13 CG iterations). 
First we should point out that the results shown in Table 14.21 are consistent with the 
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60 




Fig. 4.4. Solution u™'" on coarsest grid (no = 64) and /3 = 10 ^; about 11.5% of constraints 
are inactive. 

Table 4.2 

Iteration count for the two-level preconditioned CGS method for f} = 10~*; iteration count for 
unpreconditioned CG is in parentheses. 



SSNM iteration 


64 


128 


256 


512 


1024 


2048 


1 


2 (13) 


4(12) 


4(13) 


4(13) 


3 (12) 


3 (12) 


2 


7(13) 


6 (12) 


5 (13) 


4(13) 


4(12) 


4(12) 


3 


6 (12) 


6 (12) 


5 (13) 


4(13) 


4(12) 


4(12) 


4 


6(12) 


6(12) 


5 (13) 


4(13) 






5 


6 (12) 













theoretical estimates of Theoreni l2.2l and Remark [231 namely the number of two-grid 
preconditioned CGS iterations decreases with hj I 0. There is one notable exception: 
for the very first iterate on the coarest mesh we used as initial guess = {1, . . . , Nq}, 
so the two-grid preconditioner is the same as for the unconstrained problem, thus very 
efficient. In order to answer the question of whether the two-grid preconditioned CGS 
is more efficient than CG we adopt the point of view that in a truly large-scale context 
(three or four dimensional problems, or when applying the operators ICj requires a 
complicated sequence of operations) the cost of the solution process is largely given 
by fine-grid matrix vector multiplications (matvecs). If we accept the number of fine- 
space matvecs as a measure of efficiency we notice that the two-grid preconditioned 
CGS solves required a total of 11 CGS iterations, hence 22 matvecs at the finest level 
of 715 = 2048, compared to 36 matvecs needed by CG. As can be inferred from the 
table, the ratio (# multigrid matvecs / # CG matvecs) is decreasing with mesh-size. 
However, the above ratio would decrease much faster if da-{{Gj^)~^ , S™) would decay 
at a faster rate, as it does in the unconstrained case. Finally we should remark that 
the algorithm performed surprisingly well considering that /3 was in fact chosen quite 
small compared to what the theory suggested (/3 ~ CVh). 
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4.3. Further remarks. We return to the numerical results from Section 14.11 
which suggest that the two-grid preconditioner is closer to being of optimal order than 
Theorem 12.21 predicts, in the following sense: if restricted to subspaces of functions 
supported away from the boundary of f2J", namely spaces of the form 

Vj"jj = spanlv^P G Vf : dist(supp((^p^), 9^") > H} , 

with H > h, then the i/^^(f2)-approximation property of '^"Li is expected to be 
of almost optimal order. To see this, consider u £ ^j,H for sufficiently large H. 
Then tt^Ij^u w on d^^"^, since the size of t^jLiU decays exponentially fast away 
from supp(m), and therefore the numerical-boundary term in (|3.12p from the proof of 
Lemma l373l can be almost be dropped. Since (|3.9I) was the source of the largest term 
in the estimate p.l4p we expect that 

iProi^y^^iGT - Mf)\^y.J ^ Ch^ . 

While is is not difficult to formalize the above argument we should point out its main 
consequence: the number of generalized eigenvalues in (T(tJ"^, A^^^) that are signifi- 
cantly far from 1 is proportional to the number of nodes supported near the boundary 
of ftj^. If the dimension d is greater than 1, then this number normally increases with 
resolution, and certainly exceeds the relatively low number of iterations we found 
in Section 14.21 This indicates that the number of generalized eigenvalues which are 
0{yJTij) away from 1 dominate the computations. While the two-level preconditioner 
does not have optimal order, the above argument together with the analysis from 
Section [3] also suggest that in order to improve the presented preconditioning tech- 
nique one has to tackle the diverging action of the two operators on the nodal basis 
functions that are supported near dVl'^^. 

5. Conclusions. We have constructed two-level preconditioners for the linear 
systems arising in the semismooth Newton solution process for a control problem con- 
strained by smoothing linear operators with non-negativity inequality constraints on 
the control. The preconditioner is similar to the one constructed by Draganescu and 
Dupont in [6] for the problem without inequality control-constraints, and it maintains 
some of its qualitative behaviour: its approximation properties improve with increas- 
ing resolution. However, even though the approximation qualities of the constructed 
preconditioner are of suboptimal order, the construction and the analysis form an im- 
portant stepping stone towards finding optimal order preconditioners for the systems 
under scrutiny. The question of their practical importance and how they compare in 
efficiency with the similar preconditioners developed by Draganescu and Petra in [7] 
for systems arising in interior point methods requires a more involved study and forms 
the subject of current research. 
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